#include <iostream>
#include <iomanip>
#include <vector>
#include <array>
#include <string>
#include <sstream>
#include <numeric>
#include <cstring>
#include "Core/VecCompare.h"
#include "Core/TensorExpr.h"
#include "FiniteDiff/Multigrid.h"
#include "Embedded/EmbeddedGrid.h"
#include "Embedded/EmLevelOp.h"
#include "Embedded/EmFuncFiller.h"
#include "Embedded/EmIntergrid.h"
#include "Embedded/EMixedDer.h"
#include "Embedded/BlockRelaxation.h"
#include "../example/common/RotatedBox.h"
#include "../example/common/TestUtility.h"

#if DIM == 2

template class Tensor<Real,2>;

const int Dim = 2;
const char *bcTypes = "DDDD";
const Real theta = M_PI / 6;
Vec<Real, 3> B { 1.0, 0.0, 2.0 };
Vec<Real, 3> A;

template <class T1, class T2>
std::array<Real,3> testMixedDer(const EmbeddedGrid<2> &emGrid, const T1 &PHI, const T2 &D2PHI)
{
  EmFuncFiller<Dim> eff(emGrid);
  EmLevelOp<Dim>    elv(emGrid);
  EMixedDer         emd(emGrid, A);

  RectDomain<Dim> rd = emGrid;
  Box<Dim> gbx = rd.getGhostedBox();
  Tensor<Real, Dim> phi(gbx);
  Tensor<Real, Dim> d2phi(gbx);
  Tensor<Real, Dim> err(gbx);
  Tensor<Real, 1> bcData(emGrid.bdryNodes.size());
  eff.fillPointValues(phi, CellCentered, PHI);
  eff.fillBdryValues(bcData, "curvilinear", PHI);
  eff.transferBdryValues(phi, bcData);
  eff.fillPointValues(err, CellCentered, D2PHI);
  emd.apply(phi, d2phi);
  elv.zeroExtNodes(d2phi);
  err = err - d2phi;
  elv.zeroExtNodes(err);

  std::array<Real, 3> errnorm;
  for(int p=0; p<=2; ++p)
    errnorm[p] = elv.NRM(err, CellCentered, p);
  return errnorm;
}

template <class T1, class T2>
void doTest(const std::vector<int> &gridSize, const T1 &PHI, const T2 &D2PHI)
{
  const Real tol = 1e-12;
  using rVec = Vec<Real, Dim>;
  rVec basis[2] {
      { cos(theta), sin(theta) },
      { -sin(theta), cos(theta) }
  };
  rVec translate {sin(theta), 0.0 };
  const int Nx = gridSize.back();
  translate = translate + 1.0 / Nx;
  const int Lx = std::ceil((cos(theta) + sin(theta)) * Nx) + 2;
  RotatedBox<Dim> rb(1.0, basis);
  const int numGrid = gridSize.size();
  auto PHI_ = [&](rVec a) {
    a = a - translate;
    rVec x { cos(theta) * a[0] + sin(theta) * a[1],
             -sin(theta) * a[0] + cos(theta) * a[1] };
    return PHI(x);
  };
  auto D2PHI_ = [&](rVec a) {
    a = a - translate;
    rVec x { cos(theta) * a[0] + sin(theta) * a[1],
             -sin(theta) * a[0] + cos(theta) * a[1] };
    return D2PHI(x);
  };
  std::vector<std::array<Real, 3>> errnorm(numGrid);
  for(int n = numGrid-1; n >= 0; --n) {
    const int N = gridSize[n];
    const int L = Lx << (numGrid - n - 1);
    RectDomain<Dim> rd(Box<Dim>(0, L-1), 1.0/N, CellCentered, 2);
    EmbeddedGrid<Dim> emGrid(rd,
                             [&](rVec q, rVec &closest, rVec &normal, bool &contained, Real tol) {
                               q = q - translate;
                               rb.check(q, closest, normal, contained, tol);
                               closest = closest + translate;
                             }, false, 0.15, tol);
#ifndef NDEBUG
    if(get_dbglevel() >= 2)
      emGrid.dump(tmpos);
#endif
    std::cout << "\n\nTesting h = " << oneOver(gridSize[n]) << std::endl;
    errnorm[n] = testMixedDer(emGrid, PHI_, D2PHI_);
  }
  printConvergenceTable(&gridSize[0], errnorm);
}

int main()
{
  reset_dbglevel(1);
  A[0] = B[0] * ipow<2>(cos(theta)) - B[1] * sin(theta) * cos(theta) + B[2] * ipow<2>(sin(theta));
  A[1] = (B[0] - B[2]) * sin(2 * theta) + B[1] * cos(2 * theta);
  A[2] = B[0] * ipow<2>(sin(theta)) + B[1] * sin(theta) * cos(theta) + B[2] * ipow<2>(cos(theta));
  std::cout << "A = " << A << "\n";
  std::cout << "============================================================================\n";
  doTest({512, 256, 128, 64},
         [](Vec<Real, 2> x) {
           return sin(2 * M_PI * x[0]) * cos(2 * M_PI * x[1]);
         },
         [&](Vec<Real, 2> x) {
           Real xx = (-4*M_PI*M_PI) * sin(2*M_PI*x[0]) * cos(2*M_PI*x[1]);
           Real xy = (-4*M_PI*M_PI) * cos(2*M_PI*x[0]) * sin(2*M_PI*x[1]);
           Real yy = (-4*M_PI*M_PI) * sin(2*M_PI*x[0]) * cos(2*M_PI*x[1]);
           return B[0]*xx + B[1]*xy + B[2]*yy;
         });
  std::cout << std::endl << std::endl;
  return 0;
}

#else

int main()
{
  std::cout << "============================================================================\n";
  std::cout << "testMixedDerTrunc is skipped for 3D. " << std::endl;
  return 0;
}

#endif // DIM == 2
